
#include <cassert>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
using namespace std;

constexpr double pi = 3.1415926535897932;
struct Cd {
  static const double coe1;
  static const double coe2;
  double operator()(double A) const { return coe1 * pow(A, 0.1) + coe2 / A; }
  double g1(double A) const {
    return coe1 * 0.1 * pow(A, -0.9) - coe2 / (A * A);
  }

  double gAlpha1(double xk, double pk, double alpha) const {
    return g1(xk + alpha * pk) * pk;
  }
  double gAlpha2(double xk, double pk, double alpha) const {
    xk = xk + alpha * pk;
    return -0.09 * coe1 * pow(xk, -1.9) + 2 * coe2 / (xk * xk * xk);
  }
};

struct WolfeCon {
  double mu1{1.0e-6}, mu2{0.99};
  Cd f;
  bool operator()(double xk, double gk, double pk, double alpha, double gkp1,
                  double xkp1) {

    if (f(xkp1) > f(xk) + mu1 * alpha * pk * gk)
      return false;
    if (abs(gkp1 * pk) > mu2 * abs(gk))
      return false;

    return true;
  }
};

struct golden {
  Cd f;
  const double g = 0.5 * (sqrt(5) - 1);

  double operator()(double xk, double &a, double &b, double pk) const {
    double x1, x2, f1, f2;
    x1 = a + (1 - g) * (b - a);
    x2 = a + g * (b - a);
    f1 = f(xk + x1 * pk);
    f2 = f(xk + x2 * pk);

    // if (0.5 * (b - a) < TOL || f1 == f2)
    //   break;
    if (f1 < f2) {
      b = x2;
      x2 = x1;
      f2 = f1;
      x1 = a + (1 - g) * (b - a);
      f1 = f(xk + x1 * pk);
    } else {
      a = x1;
      x1 = x2;
      f1 = f2;
      x2 = a + g * (b - a);
      f2 = f(xk + x2 * pk);
    }
    return 0.5 * (a + b);
  }
};

// double const Cd::coe1 = 7.5209e-3;
// double const Cd::coe2 = 29.8416e-3;

double const Cd::coe1 =
    1.2 * 0.074 * 2.05 / pow((1.23 * 35 / 17.8e-6 * sqrt(11.8)), 0.2);
double const Cd::coe2 = 0.3 * 0.3 / pi / 0.96;

double GoldenSearch(const Cd &f, double a, double b,
                    const int max_step = 1000000, const double TOL = 1.0e-11);

double LineSearch(const Cd &f, double xk, const int max_step = 1000000,
                  const double TOL = 1.0e-6);

double LineNudonSearch(const Cd &f, double xk, const int max_step = 1000000,
                       const double TOL = 1.0e-6);

double LineGoldenSearch(const Cd &f, double xk, const int max_step = 1000000,
                        const double TOL = 2.0e-6);
/**
 * @brief 函数到了右侧就是非凸的了  而很多算法都是基于凸函数
 *
 * @return int
 */
int main() {
  Cd f;
  double solution = GoldenSearch(f, 1e-5, 1e2);
  cout << setprecision(15) << scientific;
  cout << "x = " << solution << "\nCd = " << f(solution) << endl;

  cout << Cd::coe1 << "\n" << Cd::coe2 << "\n";
  ofstream osf{"cd.plt"};
  for (int i = 1; i < 1000000; ++i) {
    osf << i * 0.1 << "  " << f(i * 0.1) << "\n";
  }
  osf.close();

  double solution2;
  solution2 = LineSearch(f, 40);
  cout << "linesearch x2 = " << solution2 << "\nCd2 = " << f(solution2) << endl;

  double solution3;
  solution3 = LineNudonSearch(f, 10);
  cout << "linesearch x3 = " << solution3 << "\nCd3 = " << f(solution3) << endl;

  double solution4;
  solution4 = LineGoldenSearch(f, 100000);
  cout << "linesearch x4 = " << solution4 << "\nCd4 = " << f(solution4) << endl;
  return 0;
}

double LineGoldenSearch(const Cd &f, double xk, const int max_step,
                        const double TOL)

{
  const double tol2{1.0e-14};
  double alpha{1.0}, gk, xkp1, pk, gkp1, min_alpha, max_alpha{};
  WolfeCon condition;
  double alphakp1;
  bool tag{false};
  int step{0};
  golden gg;

  do {
    ++step;
    tag = false;

    min_alpha = 1.0e-7;
    max_alpha = 1.0e6;

    gk = f.g1(xk);
    if (abs(gk) < tol2)
      break;
    if (gk > 0) {
      pk = -1.0;
      max_alpha = xk - 1.0e-1;
    } else {
      pk = 1.0;
    }

    // 黄金分割子迭代
    while (tag != true) {
      // f.galpha1 最后会越来越接近0  会产生奇异问题
      alpha = gg(xk, min_alpha, max_alpha, pk);
      xkp1 = xk + alpha * pk;
      gkp1 = f.g1(xkp1);
      if (abs(gkp1) < tol2)
        break;
      tag = condition(xk, gk, pk, alpha, gkp1, xkp1);
    //   cout << "alpha = " << alpha << endl;
    }

    xk += alpha * pk;
    cout << "alpha = " << alpha << endl;
    cout << "xk = " << xk << endl;

    //
  } while (alpha > TOL && step < max_step);
  cout << "step =" << step << endl;

  return xk;
}

double GoldenSearch(const Cd &f, double a, double b, const int max_step,
                    const double TOL) {
  const double g = 0.5 * (sqrt(5) - 1);
  double x1, x2, f1, f2;
  x1 = a + (1 - g) * (b - a);
  x2 = a + g * (b - a);
  f1 = f(x1);
  f2 = f(x2);
  for (int i = 0; i < max_step; ++i) {
    if (0.5 * (b - a) < TOL || f1 == f2)
      break;
    else if (f1 < f2) {
      b = x2;
      x2 = x1;
      f2 = f1;
      x1 = a + (1 - g) * (b - a);
      f1 = f(x1);
    } else {
      a = x1;
      x1 = x2;
      f1 = f2;
      x2 = a + g * (b - a);
      f2 = f(x2);
    }
  }
  return 0.5 * (a + b);
}

// double LineGoldenSearch(const Cd &f, double xk, const int max_step = 1000000,
//                         const double TOL = 1.0e-11) {}
// 函数不是凸的  会有问题
double LineNudonSearch(const Cd &f, double xk, const int max_step,
                       const double TOL) {
  const double tol2{1.0e-14};
  double alpha{1.0}, gk, xkp1, pk, gkp1;
  WolfeCon condition;
  double alphakp1;
  bool tag{false};
  int step{0};

  do {
    ++step;
    alpha = 1.0;
    tag = false;

    gk = f.g1(xk);
    if (abs(gk) < tol2)
      break;
    if (gk > 0)
      pk = -1.0;
    else
      pk = 1.0;

    // xkp1 = xk + alpha * pk;
    // tag = condition(xk, gk, pk, alpha, f.g1(xkp1), xkp1);
    double grad_alpha2;
    // 牛顿子迭代
    while (tag != true) {
      // f.galpha1 最后会越来越接近0  会产生奇异问题
      grad_alpha2 = f.gAlpha2(xk, pk, alpha);

      alphakp1 = alpha - f.gAlpha1(xk, pk, alpha) / grad_alpha2;
      alpha = alphakp1;
      xkp1 = xk + alpha * pk;
      gkp1 = f.g1(xkp1);

      if (abs(gkp1) < tol2)
        break;
      tag = condition(xk, gk, pk, alphakp1, f.g1(xkp1), xkp1);
    }

    xk += alpha * pk;
    cout << "alpha = " << alpha << endl;
    cout << "xk = " << xk << endl;

    //
  } while (alpha > TOL && step < max_step);
  cout << "step =" << step << endl;

  return xk;
}

// /*  最慢
double LineSearch(const Cd &f, double xk, const int max_step,
                  const double TOL) {
  const double tol2{1.0e-15};
  double alpha{1.0}, gk, xkp1, pk, gkp1, dx;
  WolfeCon condition;

  bool tag{false};
  bool decreased{false};
  bool increased{false}, gk_equal0{false};
  //   gk = f.g1(xk);
  int steps{0}, substep;
  do {
    gk = f.g1(xk);
    if (abs(gk) < tol2)
      break;

    ++steps;

    if (gk > 0)
      pk = -1.0;
    else
      pk = 1.0;

    // cout << pk << endl;
    // { return 1.0; }
    decreased = increased = tag = false;
    alpha = 1.0;
    // dx = 0.0;
    // 迭代子步  确定步长 alpha
    substep = 0;
    while (true) {
      substep++;
      //   cout << "substep = " << substep << endl;
      //   cout << "alpha = " << alpha << endl;
      //   cout << " xk = " << xk << endl;

      gk = f.g1(xk);
      xkp1 = xk + pk * alpha;
      gkp1 = f.g1(xkp1);
      tag = condition(xk, gk, pk, alpha, gkp1, xkp1);

      //   途中发现有导数为0
      if (abs(gkp1) < tol2) {
        xk = xkp1;
        gk = gkp1;
        // gk_equal0 = true;
        break;
      }
      // 满足 wolf 强条件
      if (tag == true) {
        //   之前减小过
        if (decreased)
          break;
        else {
          // 跨度太大，两者异号 但也满足 wolf条件
          if (gkp1 * gk < 0.0)
            break;

          alpha *= 2;
          increased = true;
        }
      }
      //  不满足 wolf 条件
      else {
        alpha *= 0.5;
        decreased = true;
        // 之前增加过
        if (increased)
          break;
        else
          continue;
      }
    }
    xk += pk * alpha;
    // cout << "substep = " << substep << endl;
  } while (abs(gk) > tol2 && alpha > TOL && steps < max_step);
  cout << "step =" << steps << endl;
  return xk;
}
// */